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Abstract 

We suggest a general, simplified, computationally efficient approach to simulate 
high dimensional Lindblad equations. The approach is based on a low-rank approx- 
imation of the density. The practical efficiency is demonstrated on the example of 
oscillation revivals of atoms interacting resonantly with a slightly damped coherent 
quantized field. 

1 Introduction 

Numerical simulations of high dimensional Lindblad equations are routinely performed 
using ensemble averages of quantum Monte-Carlo trajectories [2, 9, 5]. We propose here 
another approach, which consists in approximating the evolution of the n x n density 
matrix p solution to the differential Lindblad equation using a reduced dynamics on the 
set of density matrices of some fixed rank m < n. This reduced dynamics is obtained by 
taking the orthogonal projection of onto the tangent space to this set of matrices of 
rank m. 

Such an approximation strategy, based on orthogonal projections onto low dimensional 
manifolds, has already been proposed in [1, 7] in the context of quantum filtering. The 
goal then was to construct a reduced order quantum filter for a spin-spring system. The 
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submanifold on which the dynamics was projected was the (real) 4-dimensional manifold 
constructed with the tensor products of arbitrary two-level states and pure coherent states. 
In [1], a systematic approach to approximate solutions to matrix Riccati differential equa- 
tions by low-rank matrices has been likewise developed. 

We combine here the ideas of [4] and [1] to approximate arbitrary Lindblad equations. 
The proposed approximation depends only on one adjustable parameter, the low-order 
rank m. When the size n x n oi the density matrix exceeds the computing possibilities 
available -which is often the case in practice even for rather simple physically relevant 
systems-, such an approximation can be very useful to compute the (approximate) time 
evolution of p from a known initial value with rank smaller than or equal to m. This 
approximation strategy is tested here on oscillation revivals of Na atoms in a slightly 
damped mesoscopic field with n photons in average [8]. In a first stage, we consider the 
original half-spin/spring case [3] with Na = 1 and n = 15. For a photon lifetime in the 
order of one hundred vacuum Rabi periods, we compare the numerical solution p{t) to the 
Lindblad equation (computed for reference) with different numerical low-rank solutions 
from m = 2 to 6. Our results show a very satisfactory agreement. We next take the much 
larger values Na = 50 and n = 200. The density matrix is then too large for an explicit 
comparison with the direct numerical integration of the Lindblad differential equation to 
be possible using a standard computer {n ~ 15000). For a photon lifetime in the order of 
ten thousands vacuum Rabi periods, we have checked that our numerical solution using 
a reduced model of rank m = 8 is consistent with the analytic damping model proposed 
in [8]. We also observe numerically that the larger the rank m, the more accurate the 
approximation. 

Our article is articulated as follows. We derive in Section 2 the low rank dynamics 
approximating the Lindblad equation. We describe the numerical scheme we use in Sec- 
tion 3. In both sections, our general purpose strategy is specifically adjusted to account 
for the fact that, usually, the Hamiltonian part dominates the decoherence part. Section 4 
presents our numerical experiments on oscillation revivals. We draw some conclusions in 
Section 5. 

Acknowledgements The authors thank Silvere Bonnabel, Michel Brune, Michel De- 
voret, Jean-Michel Raimond for useful references and stimulating interactions. The authors 
were partially supported by the ANR, Projet Blanc EMAQS ANR-2011-BS01-017-01. 

2 The Low rank differential equations 

Consider a Lindblad equation with, for simplicity (see last paragraph of this section), a 
single decoherence operator L, 

j^p = ~z[H, p] - liL^Lp + pL^L) + LpL\ (1) 

where p is n x n non- negative Hermitian matrix with Tt (p) = 1, H is a n x n Hermitian 
matrix and L is a n x n matrix. Our purpose is to approximate, for n large, the above 
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dynamics on the set of non-negative Hermitian matrices of rank m, m being an integer 
presumably much smaller than n, prescribed beforehand. We now formalize this. 

A density p of rank m < n can be decomposed as 



P 



UaU^ (2) 
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where a is a m x m strictly positive Hermitian matrix, U a. n x m matrix with U'^U = 1 
and Im of course denotes the m x m identity matrix. The set of non-negative Hermitian 
matrices of rank m and trace one can be seen as a sub-manifold, denoted by T>m, of the 
Euclidean space of n x n Hermitian matrices equipped with the Frobenius scalar product. 
For p G Vm, given by (1) does not belong in general to the tangent space to Vm at p. 
The rank is therefore not necessarily preserved by the evolution governed by (1). To correct 
for this, the most natural option is to consider, for any p G Vm, the orthogonal projection 
of -^p given by (1) onto the tangent space to Dm at p. Denoting by H^ the projection 
operator, we therefore consider the differential equation 

|p = {-t[H, p] - \{L^Lp + pL^L) + LpLt) (3) 

set on which can be seen as a rank m approximation of (1). 

Making the approach practical requires to now give an explicit formulation of such 
a rank m approximation. A lifting procedure, inspired from [1] and described in some 
more details below, consists in introducing two coupled differential equations for U and a 
corresponding to the generic decomposition (2) for p G Vm'- 

= -lAU + (I„ - UU^) {-i{H -A)- \L^L + LUaU^ L^Ua'^U^) U, (4) 

= -i[U\H - A)U, a] - \{U^L^LUa + aU^L^LU) + U^LUaU^L^U 

+ i^Tr ((Lt(I„ - f/f/t)L UaU^) 1^. (5) 

In (4)- (5), A denotes an arbitrary Hermitian operator that may depend on time. Then 
standard manipulations, which we omit here for brevity, show that 

• U'^U remains equal to Im, 

• a remains Hermitian, positive and of trace one, 

and that the evolution of p = UaW then solves (3), which in this particular instance reads 

d 



dt 



p = -i[H, p] - \{L^Lp + pLtL) + LpLt - (I„ - P,)LpLt(I„ - P^) 



+ -<^f^^^P^, (6) 



where Pp = UW only depends on p since it corresponds to the orthogonal projection on 
the image of p. 
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Notice that (6) does not depend on the arbitrary matrix A: its entries can be seen as 
gauge degrees of freedom. The specific choice A = H yields 

= -iHU + {In - UU^) {-^L^L + LUaU^L^Ua-^U^) U, (7) 
= -^(U^L^LUa + aU^L^LU) + U^LUaU^L^U 

dt 

+ ^Tr ((Lt(I„ - UU^)L UaU^) I^, (8) 

where we note that H only appears in the dynamics for U and not in the dynamics of a, 
a choice that is particularly appropriate when H dominates L. In that case indeed, (8) 
may be understood as a slow evolution as compared to the dynamics (7). The efficiency of 
our numerical procedure, described in the next section, will significantly benefit from this 
particular decomposition. 

We now explain how we obtain (4) and (5). A nxn Hermitian matrix ^ in the tangent 
space aX p = U aU'^ to Dm admits the parameterization 

e = z[r/,p] + f/^f/^ = f/(z[f/^r/f/,(T]+^)f/t (9) 

where rj is any nxn Hermitian matrix and q any mxm Hermitian matrix with zero trace. 
This results from the definition of the tangent map of the submersion (f/, a) H- UaW with 
the infinitesimal variations 6U = irjU and 5a = The parameterization (9) is onto, but 
not one-to-one: note that different rj and q may indeed yield the same C,- The projection 
n^(J^p) corresponds to the tangent vector ^ associated to r] and ^ minimizing 

Tr (^I^Sp + pS^ + LpL^ - i[r],p] - UqU^Y^ , 

where 5* = —iH — The first order stationary conditions versus ^ and r] read 

[Sp + pS^ + L'^pL - i[r], p] - UqU^) U = \Im , 
[Sp + pS^ + L'^pL - z[r/, p] - UqU\ p] = 0, 

where the real scalar A is implicitly given by the constraint Tr [q] = 0. The first equation 
yields 

<^ = U^[Sp + pS^ + L^pL - t[v,p])U- ^Tr {U^{Sp + pS^ + pL)U) I„. (10) 

We insert this value of q into the second equation and obtain, after some easy manipulations 
using U'^U = Im, 

(I„ - P) [Sp + LpLt - n^p)p = p{pS^ + LpLt + ipri) (I„ - P) 

where P = UW is the orthogonal projector on the image of U. This matrix equation 
admits the following general solution 

7] = -i{In - P){SP + LpL^Ua-^U^) +i{PS^ + Ua-^U^LpL^){In - P) 

-PAP-{In-P)A{In-P) (11) 



4 



where A is an arbitrary Hermitian operator. Then, using WP = f/^ and PU = U, <; given 
by (10) reads 

? = -i[U^HU, a] - ^{U^L^LUa + aU^L^LU) + U^LUaU^L^U 

+ ^Tt ((Lt(I„ - f/f/t)L f/at/t) + tp^AU, a] (12) 

The derivatives -j^U = irjU and = ^ respectively give (4) and (5). 

We finally note that, for clarity, we have deliberately restricted our exposition to the 
case of a single decoherence term in (1). The generalization to an arbitrary number of 
decoherence terms 

= -^H^ P] + L^pLl - \{LlLyp + pLlLy), 
is straightforward: (7) and (8) become 

jU = -iHU + (In - UU^) ~\LlL, + L,UaU^LlUa~^U^^ U 



3 Numerical integrator 

For a positive integer k and the timestep we denote by Uk and at the numerical approx- 
imations of U{k6t) and a{k6t) solutions to (7) and (8) respectively. The evolution from 
time k6t to time {k + l)6t is split into the following three steps: 

• The first step consists in an approximation of the free Hamiltonian evolution f/fc+i = 
^-{iSt/2)Hjj^^ For the simulations of Section 4, we choose a (formal) third-order ex- 
pansion: 

UkH = (in -fH- fH' + z%H^) Uk. (13) 

We however note that, of course, other choices are possible. In particular, time inte- 
grators more adapted to the time discretization of the Schrodinger may be employed 
(see [6] and references therein). We will not proceed in this direction in the present 
work and, should need be, hope to return to this issue in future works. 

• The second step consists in updating both U and a now accounting for L; we set 

Uk+i = Uk+i + 5t{\ - Uk+iUl^) {-iL'^LUk+i + Lf/fc+.afct/t_^,Ltf/,+.a,"i) (14) 
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and, in two stages, 



CTfc+i = ak + St Ul^^iLUk+^(rkUl_^_iL^Uk+i^ 



followed by 



Tr Uul^L^LUk^ - UlMUk^Ul.LUk^)ak) Im 
+5t — ^ ^ ^ ^ , (15) 



(I- - K.^^Lt/.+i)a.+.(I„ - f /7;^LtLt/,+.) 

^^^+1 = — 7 iri iri (^^^ 

IV [(I^ - 'pl_^,VLUk^)ak^{lm - ^Ul^MLUk+^)) 

We note that this two-stage update corresponds to a slight modification of the exphcit 
Euler scheme so that it preserves both positiveness and the trace. In particular, the 
motivation for using (16) can be understood from the following considerations (which 
can be adapted to many similar contexts). Momentarily omitting its rightmost term 

for simplicity, we observe that, with obvious notation, (8) is of the form — a = 

at _ 

W a + aW'^ + Z a Z\ Introducing the auxiliary variables a = B a B'^ and Z = 

B Z B^^ with B solution to —B = —BW, we see that —a = Z b Z\ an equation 

at at 
that can be simply integrated using the first order forward Euler scheme Ok+x = 

^fc+i/2 Cfc+i/2 ^k+1/2- Expressed in the original unknown a, this scheme 

gives the numerator of (16) up to 6t^ terms and automatically preserves positiveness. 

The preservation of the trace is then ensured by the normalization in (16). 

The third step is similar to the first step, the third-order approximation of Uk+i = 
^-{iSt/2)Hjj^^^^ this time followed by an orthonormalization (a procedure formally 

denoted by T) to ensure U^j^-JJk+i = Im: 

Uk+i = T((\-fH- fH' + t%H') Uk^i). (17) 



Three comments, different in nature, are in order. 



Computational cost The above numerical scheme essentially uses right multiplications 
of H, L, L"^ hj n X m matrices, as, for example, the products HU, H^U = H{HU), LU, 
L'^{LU). To evaluate such products does not require string n x n matrices since usually 
H and L are defined as tensor products of operators of small dimensions. When n is very 
large and m is small, this point is crucial for an efficient numerical implementation of the 
approach. In particular evaluations of products like HU or LU can be easily parallelized. 
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Formal error estimator To obtain an empirical estimate of the error committed when 
using the low rank approximation, we may follow [4] and compute the Frobenius norm of 
p = ^p and p± = p — n^(p) at each time step. If the Frobenius norm of p± is much 
smaller than that of p, then the approximation is considered valid. From (6), we have 
the following general formulae for any p G Vm, any Hermitian operator H and any not 
necessarily Hermitian operator L: 

- i(LtLp + pL^L) + LpLt) = -|(LtLp + pL^L) + LpL^ 

-{i.-p,)LpL\i..-p,) + ^<'^''!:-'^'^ p, 

Thus we have, 

n = iln - P,)LpL\K - P,) - '^i^f^^l^P^ (18) 

where Pp = UW since p = UaU^ . Classical computations using standard properties of 
the trace show that, at each time step, Tr (p^) and Tr (p^) may be numerically evaluated 
with a complexity similar to the complexity of the numerical scheme. There is no need to 
explicitly compute p and pj_ as n x n matrices before taking their Frobenius norm. 

Choice of the initial condition Initial conditions for a and U need to be deduced from 
a given initial condition po. When the rank of po is larger than or equal to m, we form (Tq 
as the diagonal matrix consisting of the largest m eigenvalues of po with sum normalized 
to one and we form Uq using the associated eigenvectors. When the rank of po is strictly 
less than m then we proceed as follows. 

Assume in a first step that po is a pure state \ipo){'ipo\- It is then natural to take for (Tq 
a diagonal matrix where the first diagonal element is 1 — (m — l)e and the over ones are 
equal to e. Here e is a positive number much smaller than 1 (typically 10~^). Then Uq is 
constructed, up to an orthonormalization preserving the first column, with as the first 
column, Hlipo) as the second column, . . . , H"^~^\iIjq) as the last column. 

When the rank of po is strictly larger than 1 (but still strictly less than m), one can 
easily imagine the following mixed procedure: put the non-zero eigenvalues in the first 
diagonal elements of cxo, their associated eigenvectors in the first columns of Uq, set the 
remaining diagonal elements of ctq to e and complete the remaining columns of Uq by 
iterates of H on these eigenvectors. 

4 Numerical tests for oscillation revivals 

The collective behavior of Na atoms resonantly interacting with a quantized single-mode 
field trapped in an almost perfect cavity is modeled by the following Lindblad master 
equation (see, e.g., [8]) 

= ^ [at J- _ aJ+, p] - K{np/2 + pn/2 - apat) (19) 
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where f2o > is the couphng strength (vacuum Rabi pulsation), a the field annihilation 
operator, n = a^'a the field photon-number operator and k, > the inverse of the single 
photon life-time. The Na atoms are described here with the spin J = Na/2 representation. 
It involves the Na + I Dicke states \ J, -Na/2), \J,{-Na + l)/2), \ J, Na/2). These 
Dicke states are labelled by the index fi = to fi = Na and are denoted by for 
/i = all atoms are in the same ground state; for fi = Na + 1 all atoms are in the same 
excited state. With this notation, the atomic lowering operator, J~ = {J~^y , is given by 

J-|/i) = y/l2{Na- I2)\l2-1), fI = 0,...,Na. 

For Na = I atom initially in the excited state, a field initially in a coherent state 
with n = 15 photons (truncation to 30 photons), and a damping factor k = Qq/500, 
we have compared, until the first revival appearing around t = low-rank tra- 

jectories, solutions to (7,8), with the original full-rank trajectory solution to (1) where 
H = i^{a^J^ — aJ+) and L = y/na. In this simple case performing the full-rank simula- 
tion for the sake of comparison is indeed feasible. Initializations of U and a are performed 
according to the procedure explained at the end of Section 3. 

Figures 1, 2 and 3 show that ranks m > 4 ensure a fidelity higher than 98% over the 
interval of simulation, and that the fidelity increases with the rank. This is corroborated 
by the fact that, for m = 4 and m = 6, the Frobenius norm of p± is always less than 1% 
of the Frobenius norm of p. We recall that the fidelity between two density operators pa 
and pb is defined by Tr {^/^/mm/K) ■ 
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Figure 1: Oscillation revival for Na = 1 and a coherent initial field with n = 15 photons; 
numerical solutions of original Lindblad master equation (1) and of the rank 2 approx- 
imation (7,8) with the adimensionalized time = flot/2\/E; (a) and (b) correspond to 
the excited atomic populations; (c) and (d) show the first 2 eigenvalues; (e) is the fidelity 
between the Lindblad solution and the rank 2 solution; in (f) the solid and dashed curves 
correspond respectively to the Frobenius norms of p and p± as explained in (18). 
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Figure 2: Oscillation revival for Na = 1 and a coherent initial field with n = 15 photons; 
numerical solutions of original Lindblad master equation (1) and of the rank 4 approx- 
imation (7,8) with the adimensionalized time = flot/2\M; (a) and (b) correspond to 
the excited atomic populations; (c) and (d) show the first 4 eigenvalues; (e) is the fidelity 
between the Lindblad solution and the rank 4 solution; in (f) the solid and dashed curves 
correspond respectively to the Frobenius norms of p and p± as explained in (18). 
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Figure 3: Oscillation revival for Na = 1 and a coherent initial field with n = 15 photons; 
numerical solutions of original Lindblad master equation (1) and of the rank 6 approx- 
imation (7,8) with the adimensionalized time = flot/2\/E; (a) and (b) correspond to 
the excited atomic populations; (c) and (d) show the first 6 eigenvalues; (e) is the fidelity 
between the Lindblad solution and the rank 6 solution; in (f) the solid and dashed curves 
correspond respectively to the Frobenius norms of p and p± as explained in (18). 
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For Na = 50 atoms initially in the same excited state and a field initially in a coherent 
state with n = 200 photons a direct numerical integration of the Lindblad master equa- 
tion (19), when k > 0, is impossible on the limited computing facilities we have access to. 
The size n of any finite dimensional approximation of p by an n x n matrix necessarily 
exceeds 10000. For the purpose of validating our reduced model, we therefore proceed oth- 
erwise. With a truncation to 300 photons (ra = 51 x 301), we report here two simulations. 
On Figure 4, we show a reference simulation that is a simulation of the exact Schrodinger 
equation: the parameter k in (19) is put to zero. The complete revival is maximum. 
On Figure 5, we simulate our reduced model, based on (7,8) with rank m = 8, for the 
parameter value k = log(2)no/(47r?7,^/^). This specific choice of k > corresponds to a the- 
oretical reduction by a factor 2 of the complete revival appearing at the adimensionalized 
time = Qot/2y/E = 2ti as predicted by the formula (45) of [8]. Our simulation is per- 
formed with the numerical scheme described in Section 3 with a timestep 5t = 1/ {Vto^/fiNa)- 
Comparing Figure 5 to Figure 4, we recover numerically that the revival amplitude around 
= 27r is indeed reduced by a factor 2. This clearly validates our reduced model. Note that 
we have checked that a smaller timestep and an higher rank m = 10 do not significantly 
change the numerical results. For completeness, we mention that the large scale numerical 
results of Figures 4 and 5 come from computations performed using a simple Matlab code 
available upon request from the second author. The computations were executed on a Dell 
Precision M4400 computer equipped with Intel(R) Core(TM) duo CPU T9600 at 2.80GHz 
with 8.00 Go RAM. The rank 8 simulation of Figure 5 typically takes about 17 hours. 

5 Conclusion 

For the numerical integration of high dimensional open quantum systems, we have proposed 
an approximation approach based on an orthogonal projection onto the set of rank m den- 
sity matrices. A numerical integration scheme of the system of equations obtained for the 
reduced, low-rank model has been suggested. Its derivation has been specifically adapted 
to the situation where the Hamiltonian evolution is must faster than the decoherence evo- 
lution. Other situations could be similarly investigated. The scheme has been implemented 
and tested on oscillation revivals of atoms interacting resonantly with a slightly damped 
coherent quantized field. The results obtained show the good quality of the approximation, 
along with an evident significant speed-up with respect to the direct simulation of the full 
system. Although definite conclusions are yet to be obtained on cases of much larger sizes 
and extensive comparisons with the classically employed Monte-Carlo approaches are yet 
to be performed, the satisfactory results obtained to date show the promising nature of the 
approach. We additionally note that the strategy of approximation developed here is not 
restricted to quantum optic systems. We trust that the approach can be useful for other 
high dimensional open quantum systems. 
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Figure 4: Oscillation revival for Na = 50 and a coherent initial field with n = 200 photons; 
evolution with the adimensionalized time = flot/2^/E of the excited atomic population 
(all atoms in the same excited state); numerical solution to the Schrodinger equation (19) 
with K, = 0. 
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Figure 5: Oscillation revival for A^^ = 50 and a coherent initial field with n = 200 photons; 
numerical solution of rank m = 8 approximation (7,8); we show the excited atomic popula- 
tion (all atoms in the same excited state) versus the adimensionalized time = QQt/2y/E; 
the parameter k, = log(2)f2o/(4vrn^/^) is adjusted according to [S] in order to get a reduction 
by a factor 2 of the revival oscillations around = 27r. That reduction is indeed observed 
numerically. 
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